REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis 
Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a 
collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

28-04-2010  Journal  Article 

3.  DATES  COVERED  (From  -  To) 

4.  TITLE  AND  SUBTITLE 

Thermodynamic,  Dynamic,  and  Structural  Properties  of  Ionic  Liquids  Comprised  of  1- 
butyl-3-methylimidazolium  Cation  and  Nitrate,  Azide,  or  Dicyanamide  Anions 
(Preprint) 

5a.  CONTRACT  NUMBER 

FA8650-09-M-2036 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Dmitry  Bedrov  &  Oleg  Borodin  (Wasatch  Molelcular) 

5d.  PROJECT  NUMBER 

5f.  WORK  UNIT  NUMBER 

300509WJ 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Wasatch  Molecular,  Inc. 

2141  St.  Mary’s  Drive 

Salt  Lake  City,  UT  84108 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFRL-RZ-ED-JA-20 1 0-25 1 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/RZS 

5  Pollux  Drive 

Edwards  AFB  CA  93524-7048 

10.  SPONSOR/MONITOR’S 
ACRONYM(S) 

11.  SPONSOR/MONITOR’S 
NUMBER(S) 

AFRL-RZ-ED-JA-20 1 0-25 1 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited  (PA  #10241). 

13.  SUPPLEMENTARY  NOTES 

For  publication  in  the  Journal  of  Physical  Chemistry  B. 

14.  ABSTRACT 

Molecular  dynamics  simulations  ionic  liquids  comprised  of  [bmim]  cation  and  nitrate  [N03],  azide  [N3],  or  dicyanamide 
[N(CN)2]  anions  were  conducted  using  the  polarizable  APPLE&P  force  field.  Comparison  of  thermodynamic  properties  such  as 
densities,  enthalpies  of  vaporization,  and  ion  binding  energies  as  well  as  structural  correlations  obtained  from  simulations  at 
atmospheric  pressure  and  temperature  range  298-393  K  showed  that  IL  with  the  N(CN)2  anion  shows  significantly  different 
characteristics  compared  to  ILs  with  the  N3  and  N03  anions.  This  trend  is  further  manifested  in  dynamical  properties 
characterized  by  self-diffusion  coefficients  and  molecular  rotational  relaxation  times.  We  also  examine  the  dynamic  correlations 
between  the  ions’  translational  and  rotational  motions  as  well  as  discuss  the  anisotropy  of  the  latter. 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE 

OF  ABSTRACT 

OF  PAGES 

PERSON 

Dr.  Ghanshyam  Vaghjiani 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

SAR 

46 

19b.  TELEPHONE  NUMBER 

(include  area  code) 

N/A 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


Unclassified 


Unclassified 


Unclassified 


Thermodynamic,  dynamic,  and  structural  properties  of  ionic 
liquids  comprised  of  l-butyl-3-methylimidazolium  cation  and 
nitrate,  azide,  or  dicyanamide  anions. 


Dmitry  Bedrov  and  Oleg  Borodin 

Wasatch  Molecular  Inc.,  2141  St  Mary’s  Dr.,  Salt  Lake  City,  Utah,  84108. 

Abstract. 

Molecular  dynamics  simulations  ionic  liquids  comprised  of  [bmim]  cation  and 
nitrate  [NO3],  azide  [N3],  or  dicyanamide  [N(CN)2]  anions  were  conducted  using  the 
polarizable  APPLE&P  force  field.  Comparison  of  thermodynamic  properties  such  as 
densities,  enthalpies  of  vaporization,  and  ion  binding  energies  as  well  as  structural 
correlations  obtained  from  simulations  at  atmospheric  pressure  and  temperature 
range  298-393  K  showed  that  IL  with  the  N(CN)2  anion  shows  significantly  different 
characteristics  compared  to  ILs  with  the  N3  and  NO3  anions.  This  trend  is  further 
manifested  in  dynamical  properties  characterized  by  self-diffusion  coefficients  and 
molecular  rotational  relaxation  times.  We  also  examine  the  dynamic  correlations 
between  the  ions’  translational  and  rotational  motions  as  well  as  discuss  the 
anisotropy  of  the  latter. 
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Introduction. 


Ionic  liquids  are  extensively  investigated  for  a  number  of  applications  as  potentially 
novel  and  environmentally  friendly  materials.  Applications  of  ILs  are  extremely 
broad  including  solvents  for  synthetic  and  catalytic  applications/  lubricants/'''* 
lithium  batteries, iv-vii  actuators, vi'''x  sensors, x  reaction  media, xi  replacements  for 
conventional  solvents/''  active  pharmaceutical  ingredients/'  and  explosives  and 
propellants .x'''-xv'  One  of  the  advantages  of  ILs  is  the  possibility  to  tailor  an  IL  with 
desired  properties  through  manipulation  of  a  wide  variety  of  cation  and  anion 
chemical  structures  and  their  combinations.  While  the  empirical  knowledge  has 
been  accumulated  in  numerous  works  reported  in  the  literature,  the  fundamental 
molecular  level  understanding  of  correlations  between  the  ions’  structure  or  their 
combination  and  the  bulk  properties  of  ILs  are  still  lacking. 

In  this  work  we  employ  molecular  dynamics  (MD)  simulations  utilizing  a 
highly  accurate  polarizable  force  field  to  systematically  study  the  influence  of  the 
anion  type  on  thermophysical,  structural  and  dynamical  properties  of  l-butyl-3- 
methylimidazolium,  [bmim]-based  ILs.  Specifically,  we  focus  on  the  nitrogen 
containing  anions  such  as  nitrate  [NO3],  azide  [N3],  and  dicyanamide  [N(CN)2].  These 
ILs  have  been  explored  as  novel  energetic  materials  (explosives)  and  propellants. 
The  ILs  with  these  anions  have  been  considered  as  novel  hypergolic  materials  that 
can  ignite  spontaneously  upon  contact  with  an  oxidizer.™™  Interestingly,  it  was 
found  that  the  selection  of  anion  plays  a  crucial  role  in  determining  whether  an  IL  is 
hypergolic  or  not.  For  example,  ILs  containing  nitrate  or  azide  anions  did  not  show 
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hypergolic  behavior  upon  contact  with  common  oxidizers,  while  ILs  comprised  of 
the  same  cations  and  N(CN)2  anion  were  found  to  be  hypergolic.  The  hypergolic 
behavior  of  a  material  can,  in  principle,  depend  on  a  plethora  of  factors,  including 
kinetics  of  chemical  reactions  with  an  oxidizer,  interfacial  properties  between  the 
oxidizer  and  the  fuel,  and  processing  conditions,  etc.  All  of  these  factors  might  be 
crucial  in  defining  the  experimentally  observed  trends  in  hypergolicity  of  the  above 
mentioned  ILs.  However,  before  one  can  attempt  to  correlate  any  characteristics  of 
ILs  with  their  hypergolicity  a  fundamental  understanding  of  the  influence  of  anion 
type  on  the  basic  bulk  and  interfacial  properties  of  IL  is  needed.  Understanding  of 
the  differences  and  similarities  of  [bmim][Ns],  [bmim][N03]  and  [bmim][N(CN)2] 
bulk  properties  is  the  subject  of  this  manuscript,  while  the  interfacial  properties  of 
these  ILs  will  be  discussed  in  a  follow-up  paper.  The  selected  anions  have  similar 
sizes  (contain  only  3-5  atoms),  which  is  much  smaller  than  the  size  of  the  bmim  as 
can  be  seen  from  Figure  1  where  the  chemical  structures  of  the  ions  are  shown.  Yet, 
as  we  show  below,  despite  these  similarities,  the  ILs  containing  these  ions  have  very 
different  properties. 

Several  MD  simulations  of  some  of  these  ILs  have  been  already  reported  in 
the  literature,  xvu-xxiii  Wang  and  Voth  used  a  coarse-grained  model  to  investigate 
spatial  heterogeneities  in  alkyl  imidazolium  based  ILs  with  NO3  anion.5™'  These 
simulations  showed  that  there  is  a  significant  alkyl  tail  aggregation  in  [bmim][N03] 
liquid.  Micaelo  et  al.  developed  a  united  atom  model  for  [bmim][N03]  in  the 
framework  of  GROMOS96  force  field  and  conducted  extensive  simulations  of  this 
system  as  a  function  of  temperature.™''  These  simulations  predicted  density  and 
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viscosity  of  [bmim][N03]  that  are  in  very  good  agreement  with  experiment, 
however,  the  predicted  enthalpy  of  vaporization  of  130.2  kj/mol  at  298  K  was  about 
32  kj/mol  lower  than  the  experimental  data  later  reported  for  this  system  by 
Emel’yanenko  et  al.xxiv  Cedena  and  Maginn  used  a  fully  atomistic  model  in  their 
simulations  of  [bmim][N03].xxiii  These  simulations  also  predicted  liquid  densities  in 
very  good  agreement  with  experiment  as  well  as  provided  detailed  analysis  of  structure 
and  dynamics.  Kowsari  et  al.  investigated  transport  properties  of  [bmim] -based  ILs  with 
several  anions  including  the  [NC>3].X1X’XX  These  studies  consistently  predicted  much  slower 
dynamics  (characterized  by  self-diffusion  coefficients,  viscosities,  and  ionic 
conductivities)  compared  to  experiments  for  all  ILs  investigated.  The  only  simulations  of 
ILs  with  the  N(CN2)2  anion  that  we  are  aware  of  are  those  conducted  by  the  Steinhauser 
group  for  [emim][N(CN)2]xxv'xxvl,xxvu  in  which  the  structure,  conductivity,  and  dielectric 
relaxation  of  the  IL  were  investigated. 

Finally,  we  recently  have  reported  MD  simulations  of  variety  of  ILs  using  the 
fully  atomistic,  polarizable  APPLE&P  force  field  that  has  been  developed  and 
extensively  validated  on  a  variety  of  substances  including  ILs.xxviii  This  force  field 
showed  very  good  accuracy  in  describing  the  thermophysical  and  dynamical 
properties  for  various  ILs  using  a  transferable  set  of  non-bonded  interactions.  MD 
simulations  using  the  original  version  of  the  APPLE&P  showed  good  agreement  with 
experimental  conductivities  and  viscosities  for  [bmim][N(CN)2]  and  related  ILs, 
however,  those  simulations  predicted  density  for  [bmim][N(CN)2]  that  was  about 
2.0  %  lower  than  the  reported  experimental  values.  For  [bmim][N03],  simulations 
using  the  original  version  of  the  force  field  accurately  reproduced  reported 
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experimental  densities  as  a  function  of  temperature,  but  overestimated  the  viscosity 


by  almost  a  factor  of  two  compared  to  experiment.  To  the  best  of  our  knowledge,  no 
simulations  of  ILs  with  bmim  or  other  imidazolium-based  cations  and  N3  anion  have 
been  reported  in  the  literature. 

II.  Simulation  details. 

Subsequently  to  the  publication  of  the  original  version  of  the  APPLE&P  force  field, 
the  parameters  for  N(CN)2  and  NO3  anions  were  revised  to  improve  the  agreement 
with  available  experimental  data  in  the  liquid  phase  as  well  as  description  of  crystal 
structures  of  closely-related  ILs  containing  these  anions.5™51  In  this  work,  we  used 
the  revised  APPLE&P  potential  and  therefore  some  data  reported  here  are  slightly 
different  from  the  values  reported  in  ref  [xxviii].  Also,  recently,  the  APPLE&P  force 
field  has  been  expanded  to  include  ILs  containing  azide  anions.  While  we  are  not 
aware  of  any  experimental  measurements  on  [bmim][N3]  system,  validation  on 
closely-related  IL  ([bmmim][N3])  indicated  that  the  force  field  accurately  predicts 
the  IL  conductivity  (within  30%)  and  density  (within  1%)  compared  to 
experimental  values.  Crystal  structure  of  the  l-(2-butynyl)-3-methyl-immidazolium 
azide  was  also  found  to  be  in  a  good  agreement  with  experiments..5™1 

MD  simulations  of  bulk  [bmim][N03],  [bmim][N3],  and  [bmim][N(CN)2]  ILs 
have  been  conducted  at  298,  333,  and  393  K  and  atmospheric  pressure.  Each  system 
contained  150  ionic  pairs.  MD  simulations  were  conducted  using  the  molecular 
simulation  package  Lucretius,*5™  which  has  the  capability  to  handle  polarization 
effects.  Covalent  bond  lengths  were  constrained  using  the  velocity-Verlet  form  of  the 
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SHAKE  algorithm.™1"  The  Ewald  summation  method  was  used  for  treatment  of  long- 
range  electrostatic  forces  between  partial  atomic  charges  and  between  partial 
charges  and  induced  dipoles.  A  tapering  function  was  used  to  drive  the  induced 
dipole/induced  dipole  interactions  to  zero  at  the  cutoff  of  10  A,  with  scaling  starting 
at  9.3  A.  Induced  dipoles  were  calculated  via  a  direct  iteration  with  a  predictor 
corrector  method.  A  cutoff  of  10  A  was  used  for  all  van  der  Waals  interactions  and 
the  real  part  of  electrostatic  interactions  in  the  Ewald  summation.  A  multiple  time 
step  reversible  reference  system  propagator  algorithm5™"*  was  employed.  A  time 
step  of  0.5  fs  was  used  for  bonding,  bending,  dihedral,  and  out-of-plane  deformation 
motions,  while  a  1.5  fs  time  step  was  used  for  non-bonded  interactions  within  cutoff 
radius  of  6.0  A.  Finally,  the  non-bonded  interactions  in  the  range  between  6.0  and 
10.0  A  and  reciprocal  part  of  electrostatic  interactions  were  updated  every  3fs.  Each 
system  was  initially  equilibrated  in  the  NPT  ensemble  for  at  least  1  ns,  while 
production  runs  ranged  from  5  to  20  ns  depending  on  the  system  and  temperature. 
The  length  of  the  production  run  was  chosen  to  be  long  enough  to  allow  molecules 
to  reach  a  diffusive  regime  (when  the  ion  mean  squared  displacement  shows  a 
linear  time  dependence,  MSD[t)~  t ),  therefore  allowing  an  accurate  estimation  of 
the  self-diffusion  coefficients. 

To  calculate  the  enthalpy  of  vaporization  an  ensemble  of  150  ion  pairs  in  a 
gas  phase  has  been  simulated  using  Brownian  dynamics  simulations.  In  these 
simulations,  ions  in  a  given  (predefined)  pair  were  allowed  to  interact  only  with 
each  other,  while  any  interactions  with  other  ion  pairs  were  turned  off.  All  non- 
bonded  interactions  (van-der-Waals  and  electrostatic)  were  directly  calculated  for 
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all  possible  pair  wise  interactions  without  any  cutoff  radius.  Polarization  effects 


were  calculated  taking  into  account  only  the  electrostatic  field  created  by  the  ionic 
pair.  The  duration  of  these  simulations  was  over  Ins  allowing  accurate  sampling  of 
total  energy  per  ion  pair  in  the  gas  phase.  Finally,  to  get  the  binding  energy  of  the 
ionic  pair,  a  gas  phase  simulations  of  non-interacting  (isolated)  ions  have  been  also 
conducted  in  which  only  intramolecular  interactions  were  allowed. 

III.  Results  and  Discussion 

A.  Thermodynamic  properties.  The  [bmim][N03],  [bmim][Ns],  and 
[bmim][N(CN)2]  ILs’  density,  molar  volume  of  ionic  pair  in  the  liquid,  Vm,  enthalpy  of 
vaporization  per  ion  pair,  Hvap,  and  ion  pair  binding  energy,  E±,  are  reported  in  Table 
1  for  several  temperatures.  Examination  of  Table  I  reveals  the  following  trend  for 
densities  p[bmim][N03]>  p[bmim][N3]>  p[bmim][N(CN)2]  with  about  10%  difference  between 
the  most  and  the  least  dense  ILs.  Table  I  also  shows  that  MD  simulations  using 
APPLE&P  force  field  predict  densities  that  are  within  1%  of  the  experimental  data 
reported  in  the  literaturexxxiv'xxxv'xxxvi  (shown  in  parenthesis).  In  the  temperature 
range  investigated  the  molar  volumes  per  ionic  pair  were  found  to  have  a  linear 
temperature  dependence.  The  coefficients  of  thermal  expansion  defined  as 
a=[dV/dT)/V)  were  found  to  be  5.97xl0'4,  5.47xl0'4,  and  5.94xl0'4  for 
[bmim][N(CN)2],  [bmim][N3],  and  [bmimJfNCh],  respectively,  at  333  K. 

The  enthalpy  of  vaporization  has  been  calculated  as 

Hvap—  Egas,pair  -Euq+RT  (1) 
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where  Enq  is  the  total  energy  per  mole  of  ion  pairs  in  the  liquid,  Egas.pair  is  the  total 
ion  pair  energy  in  the  gas  phase,  R  is  the  universal  gas  constant,  and  T  is 
temperature.  We  found  that  the  enthalpies  of  vaporization  followed  a  slightly 
different  trend  than  the  densities.  While  [bmim][N(CN)2]  showed  the  lowest  HVaP, 
the  highest  Hvap  was  obtained  for  [bmim][N3].  Comparison  of  simulation  enthalpies 
of  vaporization  with  experimental  data  is  not  straight  forward  because 
experimental  enthalpy  of  vaporization  is  one  of  the  hardest  thermodynamic 
properties  to  measure  due  to  the  negligible  vapor  pressure  of  the  ILs.  Experimental 
measurements  often  have  to  be  conducted  at  elevated  temperatures  (to  attain  a 
reasonable  vapor  pressure)  where  chemical  stability  of  ILs  might  become  a  point  of 
concern.  In  Table  I  we  show  available  experimental  data  for  Hvap  for  the 
[bmim][N03]xxiv  and  [bmim][N(CN)2]xxxvii  as  obtained  using  a  combined  approach  of 
combustion  calorimetry  and  high-level  quantum  chemistry  calculations,™''  and  by 
ion  pair  fragmentation  spectroscopy.xxxviii  The  simulation  predicted  values  of  Hvap 
are  about  15  kj/mol  lower  for  both  ILs.  Similar  underestimation  of  Hvap  values 
predicted  from  simulations  using  APPLE&P  force  field  has  been  observed  in  our 
previous  studies  of  these  and  other  ILs.  However,  for  [bmim][N03],  our  simulations 
gave  a  value  for  Hvap  of  145.8  kj/mol,  which  is  in  good  agreement  with  experiment 
(162.4  kj/mol).  In  comparison,  simulations  using  the  united  atom  model  of  Micaelo 
et  al.xvii  predicted  a  value  of  130.2  kj/mol,  while  simulations  using  the  OPLS-AA 
force  fieldxviii  gave  a  value  of  117  kj/mol  when  using  the  generic  version  of  the  force 
field,  and  133  kj/mol  with  the  force  field  specifically  adjusted  for  [bmim][N03]. 
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Due  to  the  difficulty  of  measuring  the  enthalpy  of  vaporization  of  ILs,  several 
empirical  correlations  have  been  developed  to  predict  this  property.  Recently, 
Verevkin  has  examined  two  of  such  empirical  correlations  for  a  number  of  ILs.50™51 
In  that  work,  Verevkin  proposed  an  empirical  correlation  between  the  Hvap  and  the 
chemical  structure  of  IL  using  available  experimental  data  for  several  ILs.  In  his 
approach,  each  atom  type  contributes  a  certain  amount  of  enthalpy  to  the  Hvap  as 
well  as  there  are  some  additional  structural  corrections  (i.e.  for  the  ring  structure  in 
pyrrolidine  moieties  or  for  the  CF3  group).  Verevkin  used  12  different  ILs,  including 
[bmim][N(CN)2],  to  fit  these  empirical  contributions  and  then  demonstrated  that 
this  contribution  approach  works  reasonably  well  for  several  other  ILs.  The  latter 
included  [bmim][NOs]  for  which  Verevkin’s  correlation  predicted  Hvap- 169.7  kj/mol 
which  is  somewhat  higher  than  the  value  of  162.4  kj/mol  obtained  from 
combination  of  combustion  calorimetry  and  ab  initio  calculations.™''  To  the  best  of 
our  knowledge  no  experimental  data  for  Hvap  of  [bmim][Ns]  has  been  reported  in 
the  literature.  Therefore,  the  predicted  from  our  simulations  values  for  Hvap  of 
[bmim][N3]  can  only  be  compared  with  the  value  obtained  using  the  Verevkin’s 
correlation.  According  to  this  correlation  the  Hvap  for  [bmim][N3]  is  151.5  kj/mol 
which  is  5kJ/mol  lower  than  the  value  of  156.5  kj/mol  predicted  for  the 
[bmim][N(CN)2].  Quantitatively  this  value  of  Hvap  for  the  [bmim][N3]  is  very  similar 
to  the  value  of  152.5  kj/mol  predicted  from  our  simulations  at  298  K  (see  Table  I). 
However,  in  comparison  with  the  other  two  ILs,  predictions  from  Verevkin 
correlation  are  qualitatively  different  than  the  results  of  our  simulations.  In  our 
simulations,  the  [bmim][N3]  has  the  largest  value  for  Hvap,  which  is  about  22  kj/mol 
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greater  than  the  value  for  [bmim][N(CN)2].  In  contrast,  the  Verevkin  correlation 
predicts  that  Hvap  for  [bmim][N3]  should  be  the  lowest  out  of  the  three  ILs 
investigated  here  and  have  a  value  very  similar  to  that  of  [bmim][N(CN)2]. 
Implications  and  validity  of  these  qualitative  trends  are  further  discussed  below  in 
light  of  other  properties  obtained  for  these  ILs. 

Finally,  from  our  simulations  we  also  calculated  the  ion  pair  binding  energy 
defined  as 

E +—E gas, isolated  ~Egas,pair  (2) 

where  Egas, isolated  is  a  combined  energy  of  isolated,  non-interacting  with  each  other 
cation  and  anion  in  a  gas  phase.  Therefore,  while  Hvap  represents  amount  of  energy 
required  to  evaporate  an  ionic  pair  from  the  liquid  phase,  the  E+  shows  the  energy 
necessary  to  separate  (to  infinite  distance)  two  ions  in  the  gas  phase.  Table  I  shows 
that  the  E±  obtained  from  our  simulations  qualitatively  follows  the  same  trend  as 
Hvap,  i.e.  [bmim][N3]  showing  the  largest  binding  energy  and  [bmim][N(CN)2]  the 
lowest.  The  validity  of  this  trend  is  supported  by  the  ion  pair  binding  energies 
obtained  from  high-level  quantum  chemistry  (QC)  calculations  between  l-ethyl-3- 
methylimidazolium  (emim)  and  the  three  anions  investigated  in  this  work.  In  Table 
2  we  compare  binding  energies  obtained  from  QC  calculations  at  the  M05-2X/cc- 
pvTz  level  with  those  obtained  from  molecular  mechanics  (MM)  calculations  using 
the  APPLE&P  force  field.  Optimal  geometries  of  ion  pairs  from  these  calculations  are 
shown  in  the  Supplementary  Information.  The  APPLE&P  predictions  show  very 
good  agreement  with  the  QC  results  further  validating  the  accuracy  of  our  force 
field.  Also,  the  binding  energies  obtained  from  QC  calculations  for  emim-anion  pairs 
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show  the  same  trend  as  those  for  bmim-anion  pairs  reported  in  Table  1.  Taking  this 
into  account,  we  believe  that  trends  predicted  for  Hvap  using  APPLE&P  are  more 
reliable  than  those  obtained  from  Verevkin  correlation. 

B.  Structure.  To  understand  the  influence  of  anion  type  on  IL  structure,  several 
structural  correlations  were  analyzed  using  the  trajectories  obtained  from  our 
simulations.  In  Figure  2,  the  cation-cation,  cation-anion,  and  anion-anion  center-of- 
mass  radial  distribution  functions,  g(r),  are  compared  at  333  K  for  the  three  ILs 
investigated.  This  figure  clearly  shows  that  while  the  cation-cation  correlations  are 
quite  similar  for  all  three  ILs,  the  cation-anion  and  anion-anion  g[r)  for 
[bmim][N(CN)2]  are  significantly  different  from  corresponding  correlations  in 
[bmim][N3]  and  [bmim][NOs],  both  of  which  have  similar  g[r)s  to  each  other.  As  can 
be  seen  from  Figure  2b,  the  first  peak  of  cation-anion  g(r)  in  [bmim][N(CN)2]  is 
significantly  smaller  and  broader  than  the  almost  identical  peaks  in  [bmim][Ns]  and 
[bmim][N03],  indicating  a  weaker  structural  correlation  between  cation  and  anion 
centers-of-mass.  The  latter  seems  to  be  consistent  with  the  weaker  binding  energy 
and  lower  enthalpy  of  vaporization  obtained  for  [bmim][N(CN)2]  (see  Table  I  and 
discussion  above).  The  anion-anion  g(r)  for  [bmim][N(CN)2]  also  shows  noticeably 
different  behavior  compared  to  g(r)  in  [bmim][N3]  and  [bmim][N03]  as  illustrated  in 
Figure  2c.  At  short  separations  (<4.5  A),  the  anion-anion  g(r)  in  [bmim][N(CN)2]  is 
larger  than  for  the  other  two  anions  indicating  that  despite  its  larger  size  the, 
N(CN)2  anions  can  easily  approach  each  other  in  some  configurations.  Yet,  the 
position  of  the  first  peak  for  the  N(CN)2  anion-anion  g(r)  is  shifted  to  larger 
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distances  (~8  A)  compared  to  the  peak  positions  for  N3  and  NO3  anions  (~6. 5-7.0  A). 
For  [bmim][N03],  the  position  and  the  relative  heights  of  the  first  and  second  peaks 
for  all  three  g(r)  obtained  from  our  simulations  are  in  very  good  agreement  with 
those  reported  by  Cedena  and  Maginn  using  atomistic  simulations  with  non- 
polarizable  force  field.™11  Also,  our  center-of-mass  g(r)  for  [bmim][N(CN)2]  are  similar 
to  those  reported  by  Schroder  et  al.  for  [emim][N(CN)2]. 

We  continue  the  comparison  of  structural  correlations  by  examining 
atomistic  intermolecular  pair  distribution  functions.  Specifically,  in  Figure  3  we 
show  g(r)  between  hydrogens  on  imidazolium  ring  (Hi)  and  anion  atoms  that  have 
large  negative  partial  atomic  charges  and,  therefore,  experience  strong  electrostatic 
interaction  (i.e.  "hydrogen  bonding")  with  imidazolium  hydrogen  atoms.  In  Figure  1 
the  labels  and  the  values  of  partial  atomic  charges  are  given  for  atoms  for  which  g(r) 
was  calculated.  Figure  3  shows  that  Hi-Ne  g(r)  in  [bmim][N3]  and  Hi-0  g(r]  in 
[bmim][N03]  have  a  strong  first  peak  at  separation  around  2.1  A;  typical  of 
hydrogen  bonding.  The  first  peak  in  the  [bmim][N3]  Hi-N e  g(r)  is  about  30%  higher 
than  in  the  [bmim][N03]  Hi-0  g(r),  however,  this  does  not  mean  that  bmim  on 
average  makes  more  hydrogen  bonds  with  N3  than  with  NO3  anion.  IL  with  NO3  has 
a  higher  number  density  of  atoms  capable  of  forming  hydrogen  bonds  with  bmim 
(the  three  oxygens  of  NO3)  compared  to  that  in  the  [bmim][N3]  (where  only  two 
nitrogens  on  azide  ends  can  participate  in  the  hydrogen  bonding).  Note,  that  the 
partial  atomic  charge  of  the  middle  nitrogen  in  azide  has  a  relatively  large  positive 
value  and,  hence,  cannot  favorably  interact  with  bmim  hydrogens.  Therefore,  in 
order  to  compare  the  number  of  hydrogen  bonds  formed  in  each  IL,  we  calculated 
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the  coordination  number  (iVc)  within  the  first  coordination  shell  of  imidazolium 
hydrogens  defined  as  a  sphere  of  3.0  A  (corresponds  to  the  position  of  the  first 
minimum  in  Hi-0  and  Hi-Ne  g(r)  in  [bmim][NOs]  and  [bmim][N3],  respectively). 
These  calculations  reveal  that  within  a  sphere  of  3.0  A  radius,  the  imidazolium 
hydrogen  has  on  average  1.6  oxygen  atoms  of  NO3  and  1.2  nitrogen  atoms  of  N3  in 
[bmim][N03]  and  [bmim][N3],  respectively.  Although  one  needs  to  remember  that 
due  to  its  chemical  structure,  the  NO3  anion  could  easily  orient  itself  such  that  two  of 
its  oxygens  favorably  interact  with  the  same  bmim  hydrogen,  while  N3  can  only 
contribute  one  nitrogen  for  a  given  hydrogen  bond.  Finally,  the  N(CN)2  anion  has 
three  nitrogen  atoms  with  large  negative  partial  charges  as  shown  in  Figure  1  and 
hence,  in  principle,  all  of  them  can  form  hydrogen  bonds  with  bmim  .  Nevertheless, 
we  calculated  separately  the  g(r)  between  the  bmim  hydrogens  and  the  end  (Ne)  and 
middle  (Nm)  nitrogens  of  N(CN)2.  Figure  3  clearly  shows  that  while  the  Hi-Ne  g(r)  in 
[bmim][N(CN)2]  has  a  well  defined  first  peak,  similar  to  those  seen  for  N3  and  NO3 
anions,  the  Hi-Nm  g(r)  shows  no  signs  of  strong  favorable  interaction  between  anion 
middle  nitrogen  and  bmim  hydrogens.  In  this  IL,  the  bmim  Hi  on  average  has  about 
1.0  Ne  and  0.2  Nm  nitrogens  of  N(CN)2  in  a  sphere  of  3.0  A  radius. 

We  also  analyzed  the  tendency  of  alkyl  tails  of  bmim  to  segregate  into 
domains.  In  their  simulations  of  [bmim][N03]  using  coarse-grained  model,  Voth  and 
Wang  reported  a  strong  aggregation  of  end-groups  of  the  bmim  butyl  tails;  an 
indication  of  a  strong  spatial  heterogeneity  in  alkyl-imidazolium  based  ILs.xxii  We 
have  conducted  similar  analysis  and  in  Figure  4  show  g[r )  for  the  methyl  carbons  on 
the  bmim  butyl  tail  as  obtained  from  our  simulations  at  393  K.  For  the  three  ILs 
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investigated,  the  butyl  end-groups  do  not  show  any  more  structural  correlation  than 
e.g.  the  cation-anion  g[r)  shown  in  Figure  2b.  Indeed,  the  first  peak  in  the  butyl 
methyl-butyl  methyl  g[r)  reaches  values  of  1.8-1. 9  indicating  an  increased 
correlation  between  butyl  end-groups  at  short  separations,  however,  at  larger 
separations  r  >  8  A  the  g[r)  does  not  deviate  much  from  unity  indicating  a 
homogeneous  distribution  of  the  end-groups.  For  comparison  we  also  show  the 
corresponding  g{r )  extracted  from  ref.  [xxii]  which  clearly  shows  a  qualitatively 
different  behavior  despite  the  fact  that  it  has  been  obtained  from  simulations  at  a 
much  higher  temperature  (700K)  where  spatial  correlations  are  expected  to  be 
weaker  than  at  393  K.  Taking  into  account  such  qualitative  difference  between  our 
g[r)  and  those  from  ref.  [xxii],  we  believe  that  the  coarse-grained  model  used  by 
Voth  and  Wang  greatly  overestimated  the  aggregation  between  the  alkyl  tails  and, 
hence,  the  structural  heterogeneities  reported  in  that  work  are  over  exaggerated. 

We  conclude  our  analysis  of  structural  correlations  with  3-dimensional  (3D) 
density  distributions  of  anion  atoms  around  the  cation.  In  this  analysis,  the  three 
carbon  atoms  on  imidazolium  ring  defined  one  of  the  coordinate  system  planes.  The 
simulation  box  was  divided  into  lattice  cells  of  size  0.4  A  and  the  relative  density 
distribution  pi/<p>,  where  pn  is  the  local  number  density  of  atoms  of  interest  in  the  /- 
the  lattice  cell  and  <p>  is  the  average  number  density  of  these  atoms  in  the  system, 
was  calculated  by  averaging  over  all  bmim  molecules  and  the  entire  trajectory. 
First,  we  compare  these  3D  distributions  for  anion  atoms  that  form  hydrogen 
bonding  with  bmim  Hi  atoms,  i.e.  Ne  in  [bmim][Ns]  and  [bmim][N(CN)2]  and  0  in 
[bmim][N03].  At  333  K,  the  maximum  values  of  pi/<p>  were  found  to  be  25.4  and 
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40.0  for  Ne  in  [bmim][N3]  and  [bmim][N(CN)2],  respectively,  and  24.8  for  0  in 
[bmim][N03].  In  Figure  5  isosurfaces  for  pi/<p>=\ 0.0  (Fig  5a)  and  18.0  (Fig.5b)  are 
shown  for  these  distributions.  These  isosurfaces  have  relatively  high  values  of  local 
density  (compared  to  the  bulk  average)  and  therefore  represent  relatively  low  free 
energy  configurations  that  frequently  occur  during  simulations.  At  pi/<p>= 10.0 
distributions  for  all  three  ILs  are  very  similar  showing  three  well  defined  regions  of 
preferable  location  of  hydrogen  bonding  atoms,  which  are  found  to  be  in  the  vicinity 
of  bmim  Hi.  These  locations  are  also  similar  to  the  preferred  locations  of  Ntf2  anion 
oxygen  atomsxl  further  indicating  the  similarity  of  the  bmim  cation  coordination  by 
various  anions.  The  similarity  of  distributions  observed  in  Figure  5a  indicates  that 
for  a  given  cation  (bmim),  the  anion  type  does  not  influence  the  distribution  of 
preferable  locations  of  anion  atoms  that  participate  in  hydrogen  bonding.  On  the 
other  hand,  the  maximum  values  of  pi/<p>  reported  above  indicate  that  despite  the 
fact  that  <p>  of  N(CN)2  Ne  atoms  is  about  10%  lower  than  that  for  Ne  of  N3  (due  to 
difference  in  overall  density  of  [bmim][N(CN)2]  and  [bmim][N3]),  the  maximum 
value  of  pi/<p>  for  the  latter  is  almost  two  times  larger.  This  indicates  that  there  are 
certain  locations  around  the  bmim  cation  in  which  nitrogens  at  N3  ends  have 
significantly  lower  free  energy  compared  to  those  on  N(CN)2  ends.  Isosurfaces  for 
higher  density  [pi/<p>= 18.0,  Figure  5b)  show  that  the  region  above  the  bmim  N-C-N 
bend  is  less  favorable  than  the  other  two  regions  on  the  other  side  of  the  bmim  ring 
where  anion  atoms  can  form  hydrogen  bonds  without  steric  interference  from  the 
methyl  and  butyl  groups  and  where  the  maximum  values  in  pi/<p>  were  found  for  all 
anions  investigated  here. 
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3D-distributions  can  also  be  helpful  in  understanding  the  relative  orientation 
of  anions  in  the  first  coordination  shell  of  bmim.  While  for  the  relatively  symmetric 
NO3  anion  its  relative  orientation  is  not  particularly  interesting,  the  orientation  of  N3 
and  N(CN)2  anions  can  be  gleaned  from  the  analysis  of  Figure  6  where  the  3D- 
distribution  of  end  and  middle  nitrogens  are  compared  for  the  same  anion.  Figure 
6a  shows  that  the  distribution  of  Ne  and  Nm  nitrogens  of  N3  are  basically  parallel  to 
each  other  indicating  that  the  longest  dimension  of  the  N3  anion  is  almost  aligned 
along  the  Hi-Ne  hydrogen  bond.  On  the  other  hand,  isosurfaces  for  Nm  and  Ne  of 
N(CN)2  are  shifted  relative  to  each  other  as  can  be  seen  from  Figure  6b.  The  relative 
location  of  these  isosurfaces  is  consistent  with  configurations  in  which  the  N(CN)2 
anion  can  form  two  hydrogen  bonds  with  the  same  bmim;  one  with  Hi  on  the  N-C-N 
bend  and  another  with  Hi  on  the  opposite  side  of  imidazolium  ring. 

C.  Dynamical  properties.  We  begin  our  analysis  of  dynamical  properties  with 
comparison  of  self-diffusion  coefficients  ( D )  of  ions  as  a  function  of  temperature.  To 
obtain  the  self-diffusion  coefficient,  mean-squared  displacements  (MSD)  for  each 
ion  type  were  calculated  as  a  function  of  time  and  are  shown  in  Figure  7  at  393  K. 
For  all  three  ILs  at  any  given  time  t  the  anion  has  a  larger  MSD  compared  to  bmim, 
which  is  expected  taking  into  account  that  all  anions  investigated  here  are  smaller 
than  bmim  cation.  Figure  7  clearly  shows  that  anion  type  has  significant  influence  on 
the  mobility  of  ions  in  ILs.  While  ion  MSD  in  [bmim][N3]  and  [bmim][N03]  are 
similar,  with  somewhat  systematically  higher  MSD  for  the  latter,  the  ion  mobility  in 
[bmim][N(CN)2]  is  significantly  larger  than  in  the  other  two  ILs.  To  obtain  the  self- 
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diffusion  coefficients,  the  MSD(t)  were  fitted  as  a  linear  function  of  time  (t)  for  data 
where  MSD(t)  >  10.0  A2 .  The  latter  condition  ensures  that  the  fitting  is  done  for  the 
timescales  on  which  ion  motion  is  diffusive.  In  Table  3  the  ion  self-diffusion 
coefficients  obtained  in  this  way  are  given  for  three  ILs  at  three  temperatures 
investigated.  Consistent  with  MSD  shown  in  Figure  7,  the  obtained  self-diffusion 
coefficients  in  [bmim][N(CN)2]  are  significantly  higher  than  for  the  other  two  ILs.  At 
298  K,  the  mobility  of  both  ions  in  [bmim][Ns]  is  an  order  of  magnitude  slower  than 
that  in  the  [bmim][N(CN)2].  At  393  K,  the  difference  between  ion  self-diffusion 
coefficients  in  these  ILs  reduces  but  still  is  about  a  factor  of  three  different.  Ratios  of 
the  anion  and  cation  self-diffusion  coefficients  in  the  same  IL  are  also  given  in  Table 
3.  These  show  a  strong  correlation  between  the  cation  and  anion  mobilities.  The 
self-diffusion  coefficient  of  anions  is  only  15-30%  higher  than  that  for  bmim  in 
[bmim][N03]  and  [bmim][N3].  In  [bmim][N(CN)2],  the  difference  between  ions 
mobilities  is  about  40-50%.  Interestingly,  both  simulations  using  non-polarizable 
force  fields™'’™**  predicted  that  the  self-diffusion  of  NO3  in  [bmim][N03]  is  slightly 
(10-20%)  slower  than  that  of  bmim,  which  seems  surprising  taking  into  account  the 
difference  in  bmim  and  NO3  sizes. 

In  our  recent  paperxl  we  showed  that  empirical  correlations  between  the 
ions’  average  self-diffusion  coefficient  and  thermodynamic  properties  of  IL  work 
reasonably  well  for  a  wide  variety  of  ILs.  Specifically,  correlations  relating  D  with 
the  Hvap,  E± ,  and  Vm,  have  been  suggested.  In  Figure  8  we  demonstrate  one  of  these 
correlations  -\og[D)~Hvap+0.18E±  for  several  ILs  that  contain  bmim  cation  and  have 
been  simulated  using  the  APPLE&P  force  field,  including  those  simulated  in  this 
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work.  As  can  bee  seen  from  Figure  8,  the  data  for  [bmim][Ns]  follow  this  empirical 


correlation  very  well  further  indicating  that  our  relatively  high  values  of  Hvap  are 
consistent  with  the  slow  dynamics  in  this  system.  No  direct  measurements  of  the 
ion  mobility  in  [bmim][Ns]  has  been  reported  in  the  literature,  however,  there  are  a 
couple  indirect  case  of  evidence  that  the  noticeably  slower  dynamics  in  this  IL  can 
be  expected  compared  to  ILs  with  other  anions.  In  ref.  [xv],  a  value  of  404  cP  at 
298K  has  been  reported  for  viscosity  of  a  supercooled  [bmim][Ns]  liquid,  which  is 
about  a  factor  of  two  larger  than  the  experimental  data  for  [bmim][NOs].xxxv  Also,  a 
recent  study  of  [bmmim] -based  ILs  and  mixtures  with  N3  and  BF4  anionsxli  showed 
that  at  313  K,  IL  with  the  N3  anion  had  about  60%  higher  viscosity  than  the  IL  with 
the  BF4  anion.  Taking  into  account  the  consistency  of  simulation  data  with  empirical 
correlation  in  Figure  8  as  well  as  experimental  data  on  related  systems  mentioned 
above,  we  believe  that  both  thermodynamic  and  dynamic  properties  predicted  from 
our  simulations  for  [bmim][N3]  are  quite  reasonable. 

We  also  analyzed  the  rotational  dynamics  of  ions  as  a  function  of 
temperature.  In  this  analysis,  for  all  ions,  except  the  N3,  a  local  coordinate  system 
was  defined  and  a  rotational  autocorrelation  function  (ACF)  for  each  unit  vector  e,  ( 
i-{x,y,z })  defining  this  local  coordinate  system  was  calculated  as: 

^c-^(r)  -  ([*.(»)■  *,(»)])  (3) 

where  e,(0)  and  e,(t)  are  the  values  of  the  unit  vectors  at  time  zero  and  t, 
respectively,  while  brackets  denote  an  ensemble  average  over  all  molecules  of  the 
same  type  and  time  origins.  The  location  of  the  origin  and  orientation  of  the 
selected  local  coordinate  system  for  each  molecule  are  given  in  Figure  1.  For  the  N3 
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anion,  a  vector  between  two  nitrogen  atoms  on  the  molecule  ends  was  considered 
instead  of  rotation  of  the  local  coordinate  system.  The  obtained  ACFs  were  fitted 
with  Kohlrausch-Williams-Watts  (KWW)  function  given  by: 

(4) 

where  tr  is  a  relaxation  time  parameter,  parameter  /?  determines  the  degree  of 
stretching  and  characterizes  the  broadness  of  the  relaxation  process,  and  prefactor 
A  allows  us  to  account  for  the  decay  which  occurs  on  time  scales  faster  than  1  ps_1. 
Rotational  relaxation  times,  (z),  were  obtained  by  integrating  eq.  4  over  time  from 
zero  to  infinity  and  are  given  in  Table  4. 

Examination  of  Table  4  reveals  several  interesting  observations.  Rotation  of 
bmim  in  all  ILs  is  quite  anisotropic.  Relaxation  of  the  x-axis,  which  is  aligned  along 
the  largest  dimension  of  bmim  molecule,  is  about  a  factor  of  3-4  slower  than  the 
relaxation  of  the  other  two  axes  for  which  the  reorientation  of  the  butyl  tail  is  not 
required.  As  with  the  self-diffusion  coefficient,  the  slowest  rotation  of  all  axes  is 
observed  for  bmim  in  [bmim][Ns]  while  the  fastest  is  seen  in  [bmim][N(CN)2].  The 
rotation  of  the  N(CN)2  anion  is  also  quite  anisotropic  with  relaxation  of  the  x-axis 
being  about  an  order  of  magnitude  slower  than  that  for  the  other  two  axes.  On  the 
other  hand,  the  relatively  symmetric  NO3  anion  shows  very  isotropic  rotation.  As  we 
discussed  above,  the  self-diffusion  coefficients  of  cation  and  anion  in  the  ILs 
investigated  were  quite  similar  (deviating  no  more  than  50%  from  each  other). 
However,  a  strong  decoupling  of  the  cation  and  anion  rotational  dynamics  can  be 
seen  from  the  data  reported  in  Table  4.  In  [bmim][N(CN)2],  the  anion’s  slowest 
relaxation  time  (xx)  is  about  4  times  faster  than  the  slowest  rotation  relaxation  time 
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for  bmim,  while  relaxations  for  the  other  two  axes  are  different  almost  by  a  factor  of 
ten.  Similarly,  about  two  orders  of  magnitude  difference  between  the  rotational 
relaxations  of  the  cation  and  anions  is  observed  in  [bmim][Ns]  and  about  three 
orders  in  [bmim][NOs].  The  latter  is  due  to  exceptionally  fast  rotational  dynamics  of 
NO3  compared  to  other  anions  investigated  here  and  is  likely  a  consequence  of  the 
symmetric  structure  of  NO3  anion.  Cadena  and  Maginn  also  analyzed  the  rotational 
dynamics  of  ions  in  their  simulations  of  [bmimJfNOs].5™  Although  they  used  a 
somewhat  different  definition  of  the  molecular  vector  for  which  the  rotational 
relaxation  has  been  determined,  the  rotational  relaxations  obtained  in  their  work 
for  bmim  are  consistent  with  our  data  in  Table  4.  They  found  that  the  rotational 
relaxation  of  bmim  changes  from  2.5ns  at  298  K  to  about  100  ps  at  393  K  (the  latter 
value  was  interpolated  using  the  reported  values  at  363  and  423K),  while  the 
average  (over  rx,  ry,  and  rz)  rotational  relaxation  from  our  simulations  changed 
from  2.8  ns  to  120  ps  in  the  same  temperature  range.  However,  for  the  NO3  anion 
Cadena  and  Maginn  reported  rotational  relaxation  times  of  91  ps  at  298  K  and  about 
3ps  at  393  K,  while  our  simulations  predict  6  and  1.8  ps  relaxation  times  at  298K 
and  363  K,  respectively. 

It  is  also  interesting  to  compare  the  temperature  dependence  of  the 
rotational  and  translational  motion,  and  thus  to  investigate  the  coupling  of 
rotational  and  translational  dynamics  as  a  function  of  temperature.  In  Figure  9  we 
show  the  temperature  dependence  of  nD,  where  n  are  the  rotational  relaxation 
times  from  Table  4  (i=  x,  y,  or  z)  and  D  is  the  ion  self-diffusion  coefficient  from  Table 
3.  If  the  rotational  dynamics  as  a  function  of  temperature  were  to  change  the  same 
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way  as  the  translational  dynamics  (i.e.  the  activation  energies  for  these  types  of 
motion  would  be  the  same),  the  product  of  the  self-diffusion  coefficient  and 
rotational  relaxation  time  would  be  constant.  As  can  be  seen  in  Figure  9a,  this  is 
mostly  the  case  for  bmim  in  the  three  ILs  investigated.  However,  for  anions,  whose 
nD  are  shown  in  Figure  9b,  only  N(CN)2  shows  this  temperature  independence  (or 
weak  dependence)  of  nD,  while  N3  and  NO3  show  strong  temperature  dependences, 
indicating  that  these  anions  have  quite  different  temperature  scaling  for  rotational 
and  translational  dynamics.  Indeed  values  in  Tables  3  and  4  show  that  e.g.  the  self¬ 
diffusion  coefficient  of  NO3  increases  by  a  factor  of  30  upon  heating  from  298  to  393 
K,  while  the  rotational  dynamics  speeds  up  only  by  a  factor  of  3-4  in  the  same 
temperature  range.  A  similar  rotational-translational  decoupling,  albeit  to  a  smaller 
extent,  has  been  observed  in  our  simulations  of  ILs  with  another  small  anion 
bis(fluorosulfonyl)imide  (FSI),  while  a  larger  bis(trifluoromethylsulfonyl)imide 
(Ntf2)  anion  showed  good  coupling  between  rotational  and  translational  motions  in 
the  temperature  range  298  K  -393  K.xlii 

Another  interesting  information  that  can  be  obtained  from  Figure  9  is  the 
average  value  of  the  ion  MSD  on  the  time  scale  corresponding  to  the  rotational 
relaxation  time  for  a  given  axis  (MSDT=6  nD).  Figure  9a  reveals  that  on  the  bmim  zx 
time  scale,  the  cation  center-of-mass  MSDT  is  about  18  A2  indicating  that  bmim  has 
to  move  out  of  its  intermolecular  "cage”  before  rotational  relaxation  of  the  x-axis 
(see  Fig  1  for  definition)  occurs.  On  the  other  hand,  for  bmim  zy  and  rz,  the 
corresponding  MSDT^6  A2,  which  is  noticeably  smaller  than  the  dimensions  of  bmim 
molecule.  This  indicates  that  the  rotational  relaxation  of  these  axes  can  occur 
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without  significant  translational  motion  of  the  molecule,  i.e.  without  leaving  the 
cage.  Similar  analyses  for  the  anions  show  that  only  rotation  of  N(CN)2  around  its 
largest  dimension  (x-direction)  requires  any  noticeable  displacement  of  the  ion 
(MSDt»6  A2).  For  all  other  anions  and  rotational  axes,  the  rotational  relaxation  can 
occur  with  minimal  translational  motion  of  the  ion.  This  is  particularly  apparent  at 
298  K  where  MSDT  less  than  0.1  A2  are  required  for  rotational  relaxation  of  N3  and 
NO3  anions.  While  this  result  can  be  expected  for  NO3  due  to  its  symmetric  structure, 
the  result  is  somewhat  surprising  for  the  linear  N3  anion.  One  might  expect  that 
formation  of  two  relatively  strong  hydrogen  bonds  by  the  N3  ends  with  the 
surrounding  bmim  molecules  would  interfere  with  the  N3  rotation  and,  therefore, 
the  rotation  of  this  anion  would  be  coupled  to  structural  rearrangement  n  the 
system  (i.e.  translational  motion  of  ions).  This  is  clearly  not  the  case,  indicating  that 
the  formed  hydrogens  bonds  between  N3  and  bmim,  while  strong,  are  not  long  lived. 
The  latter  might  be  due  to  availability  of  other  hydrogen  bonding  sites  on  bmim  ring 
that  can  exchange  the  hydrogen  bonding  with  the  N3  nitrogen  without  relaxation  of 
intermolecular  cage. 

Finally,  the  ionic  conductivity  (X)  was  calculated  for  each  IL  using  the 
approach  described  in  ref.[xxviii]  and  is  also  given  in  Table  I.  In  the  ILs  investigated 
the  degree  of  ion  dynamic  dissociation,  as  defined  in  ref.[xxviii],  was  around  0.6  and 
therefore  trends  in  ionic  conductivities  obtained  from  simulations  strongly  correlate 
with  those  observed  above  for  the  self-diffusion  coefficients.  Table  1  also  reports 
experimental  conductivities  available  in  the  literature  for  [bmim]  [N (CN)2].xliii  A  very 
good  agreement  between  the  simulation  data  and  experiment  for  this  IL  further 
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validates  the  quality  of  the  APPLE&P  force  field  and  provides  confidence  in 


simulation  predictions  reported  here  for  which  no  experimental  data  are  available. 

IV.  Conclusions. 

Molecular  dynamics  simulations  using  the  polarizable  APPLE&P  force  field  were 
conducted  to  investigate  the  influence  of  anion  type  on  properties  of  [bmim] -based 
ILs.  We  have  focused  on  ILs  with  relatively  small,  nitrogen  containing  anions  such  as 
nitrate  [NO3],  azide  [N3],  and  dicyanamide  [N(CN)2].  We  found  that  the  IL  containing 
the  N(CN)2  anion  shows  noticeably  different  density,  heat  of  vaporization,  and  ion 
binding  energy  compared  to  the  ILs  with  N3  and  NO3  anions.  Our  simulations 
showed  that  [bmim][N3]  has  the  largest  heat  of  vaporization  while  [bmim]N(CN)2] 
has  the  smallest;  a  trend  which  is  qualitatively  different  from  the  predictions  of 
some  empirical  correlations  reported  in  the  literature,  but  consistent  with  the  high- 
level  quantum  chemistry  calculations.  We  also  found  that  the  cation-anion  center-of- 
mass  pair  distribution  function  in  the  IL  with  the  N(CN)2  showed  weaker 
correlations  than  the  other  two  ILs.  On  the  other  hand,  analysis  of  the  three- 
dimensional  density  distributions  of  anion  atoms  around  the  cation  showed  that  the 
preferable  anion-cation  coordination  is  independent  of  the  anion  type.  Analysis  of 
the  butyl  end-group  distribution  showed  a  small  tendency  for  its  aggregation 
independent  of  anion  type.  However,  the  extent  of  such  aggregation  observed  in  our 
simulations  was  found  to  be  significantly  smaller  than  what  has  been  previously 
reported  in  the  literature  for  [bmim][N03]  using  coarse-grained  simulations. 
Analysis  of  ion  dynamics  showed  that  the  highest  ion  mobility  is  observed  in 
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[bmim][N(CN)2];  consistent  with  the  lowest  heat  of  vaporization  predicted  for  this 
IL.  Rotational  dynamics  of  the  anions  was  found  to  be  decoupled  from  ions 
translational  motion,  therefore  allowing  significant  rotational  relaxation  with 
minimal  translational  motion  of  the  ions. 
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Table  1.  Density,  molar  volume  per  ion  pair,  enthalpy  of  vaporization  and  binding 


energy  of  ionic  pair  as  obtained  from  MD  simulations  using  the  APPLE&P  force  field. 
Values  in  parenthesis  are  experimental  data  available  in  the  literature. 


IL 

T 

(K) 

P 

(kg/m3) 

vm 

(cm3/mol) 

Hvap 

(kJ/mol) 

E± 

(kJ/mol) 

[bmim][N(CN)2] 

298 

1049.7  (1058a) 

195.3 

132.6  (157. 2e) 

(153. 4e) 
(174f) 

366.5 

333 

1028.2  (1041a) 

199.4 

128.9 

365.3 

393 

992.2 

206.6 

122.9 

363.2 

[bmim][N3] 

298 

1073.2 

168.7 

152.5 

389.9 

333 

1053.4 

171.8 

148.7 

388.7 

393 

1019.1 

177.6 

142.2 

386.8 

[bmim][N03] 

298 

1159.9  (1154b) 

173.3 

145.8  (162.4d) 

387.0 

333 

1135.6  (1131b) 
(1136c) 

177.0 

142.2 

386.2 

393 

1097.4  (1092b) 

183.3 

136.2 

384.6 

a)  ref.  [34]  b)  ref.  [35]  c)  ref.  [36]  d)  ref.  [24]  e)  ref.  [37]  f)  ref.  [xxxviii] 
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Table  2.  Cation-anion  binding  energies  (in  kj/mol)  as  obtained  from  quantum 


chemistry  (QC)  calculations  (at  the  M05-2X/cc-pvTz  level)  and  molecular  mechanics 
using  APPLE&P  force  field.  Binding  energies  defined  as  E isolated  ions- Epair. 


Ion  pair 

QC 

APPLE&P 

[emim][DCA],geomla 

453 

470 

[emim][DCA],geom2 

444 

466 

[emim][N3] 

508 

498 

[emim][N03] 

499 

491 

a  Two  geometries  with  similar  energies  were  investigated  for  [emim][N(CN)2]. 
Pictures  and  XYZ  coordinates  of  optimized  geometries  are  given  in  the  Supporting 
Information. 
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Table  3.  Ion  self-diffusion  coefficients  obtained  from  MD  simulations  at  various 


temperatures.  Also  shown  are  the  ratios  between  the  anion  and  cation  self-diffusion 
coefficients  and  ionic  conductivity.  Where  available  experimental  data  are  shown  in 
parenthesis. 


IL 

T  (K) 

D  (1010 

cation 

m2/s) 

anion 

Dan/Dcat 

X 

mS/cm 

[bmim][N(CN)2] 

298 

0.30 

0.41 

1.37 

8.9  (lla) 

333 

1.07 

1.50 

1.40 

28.4  (24. la) 

393 

3.21 

4.95 

1.54 

73.7 

[bmim][N03] 

298 

0.058 

0.066 

1.14 

1.7 

333 

0.29 

0.37 

1.27 

8.0 

393 

1.50 

1.90 

1.27 

32.8 

[bmim][N3] 

298 

0.020 

0.023 

1.15 

0.65 

333 

0.15 

0.17 

1.13 

4.1 

393 

1.11 

1.38 

1.24 

26.8 

a)  ref  [42] 
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Table  4.  Rotational  relaxation  times  in  ps  obtained  from  MD  simulations  at  various 


temperatures.  See  text  and  Figure  1  for  definition  of  the  x,y,  and  z  axes. 


IL 

T(K) 

TX 

cation 

Ty 

TZ 

TX 

anion 

Ty 

Tz 

[bmim][N(CN)2] 

298 

1036 

315 

332 

251 

37.3 

28.9 

333 

299 

96.6 

101 

77.4 

7.8 

6.2 

393 

90.8 

28.3 

29.1 

24.2 

1.9 

1.7 

[bmim][N03] 

298 

4588 

1808 

1925 

5.1 

5.1 

8.0 

333 

1134 

402 

387 

3.1 

3.1 

3.9 

393 

200 

75.6 

77.2 

1.7 

1.7 

1.9 

[bmim][N3] 

298 

26489 

4120 

4130 

99.9 

333 

2415 

684 

753 

26.8 

393 

243 

92.0 

91.4 

6.7 
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Figure  Caption. 


Figure  1.  Molecular  structure  of  ions  used  in  simulations.  Also  shown  are  the  labels 
and  partial  atomic  charges  for  the  selected  atoms  as  well  as  the  orientation  and 
position  of  the  local  coordinate  system  for  rotational  dynamics  analysis.  Only  x  and 
y  axes  are  shown,  the  z-axis  is  perpendicular  to  the  xy-plane. 

Figure  2.  Ion-ion  center-of-mass  radial  distribution  functions  as  obtained  from  MD 
simulations  of  three  ILs  at  333  K. 

Figure  3.  Radial  distribution  functions  between  hydrogens  on  imidazolium  ring  and 
anion  atoms  with  large  negative  partial  atomic  charges  (see  Figure  1  for  label 
notation). 

Figure  4.  Radial  distribution  functions  between  the  butyl  methyl  groups  as 
obtained  from  our  simulations  of  [bmim][N03],[bmim][N3],  and  [bmim][N(CN)2]  at 
393  K.  Also  shown  is  the  corresponding  g[r)  obtained  from  the  coarse-grained 
simulations  of  [bmim][N03]  by  Wang  and  Voth  (extracted  from  ref.  [xxii]). 

Figure  5.  3D  distributions  of  pi/<p>  isosurfaces  for  Ne  (in  N3  and  N(CN)i)  and  O  (in 
N03)  atoms  around  bmim:  a)  pt/<p>= 10.0,  b)  /?;/<ya>=18.0.  Red:  [bmim][N(CN)2], 
green:  [bmim][N3],  and  blue:  [bmim][N03]. 
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Figure  6.  3D  distributions  of  pi/<p>  isosurfaces  for  Ne  and  Nm  atoms  around  bmim:  a) 
in  [bmim]  [N3]  solid  pNe/<pNe>= 17-0  (max(pNe/<pNe>)  =  40.0),  wireframe: 
pNm/<PNm>=\0-0  (max(pNm/<PNm>)=24.0);  b)  in  [bmim][N(CN)2]  solid: 
PNe/<PNe>= 10.0  (max(pNe/<PNe>)  =  25.2),  wireframe:  pNm/<pNm>= 4.5  (max(pNm/<PNm 
>)=9.\). 

Figure  7.  Mean  squared  displacements  (MSD)  of  ions  as  a  function  of  time  obtained 
from  simulations  of  [bmim][N03],  [bmim][N3],  and  [bmim][N(CN)2]  at  393  K. 

Figure  8.  Correlation  between  the  average  ion  self-diffusion  coefficient  and  the  enthalpy 
of  vaporization  (Hvap)  and  the  cation-anion  binding  energy  E+  as  obtained  from  MD 
simulations  using  APPLE&P  force  field  for  various  ILs  with  bmim  cation  at  298  K.  Data 
include  three  ILs  investigated  here  and  those  reported  in  ref.  [xl]. 

Figure  9.  Temperature  dependence  of  the  product  of  rotational  relaxation  time  (r,-)  and 
the  ion  self-diffusion  coefficient  (D)  for  a)  cation  and  b)  anion. 
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Figure  8. 
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Supporting  Information 


Optimal  low  energy  geometries  as  obtained  from  quantum  chemistry  calculations 


for  ionic  pairs.  Several  key  atom-atom  separations  are  also  provided  (normal  font- 


QC,  italic-  APPLE&P). 


[bmim][N(CN)2]  geometry  1 
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[bmim][N(CN)2]  geometry  2 
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